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Abstract. 

The dynamic scaling of mesoscopically thick films (up to 10"^ atomic layers) 
grown with the Clarke-Vvedensky model is investigated numerically for broad 
ranges of values of the diffusion-to-deposition ratio R and lateral neighbor 
detachment probability e, but with no barrier at step edges. The global roughness 
scales with the film thickness t &s W ^ j {e + a)]^ where (3 ^ 0.2 is 

the growth exponent consistent with Villain-Lai-Das Sarma (VLDS) scaling and 
a = 0.025. This general dependence on R and e is inferred from renormalization 
studies and shows a remarkable effect of the former but a small effect of the latter, 
for e < 0.1. For R > 10“^, very smooth surfaces are always produced. The local 
roughness shows apparent anomalous scaling for very low temperatures {R < 10^), 
which is a consequence of large scaling corrections to asymptotic normal scaling. 
The scaling variable R^^^ (e-|-a) also represents the temperature effects in the 
scaling of the correlation length and appears in the dynamic scaling relation of 
the local roughness, which gives dynamic exponent z 3.3 also consistent with 
the VLDS class. 


PACS numbers: 68.55.-a, 68.35.Ct, 81.15.Aa , 05.40.-a 
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1. Introduction 

Stochastic modeling of molecular-beam epitaxy (MBE) attracted much interest in the 
last decades due to the importance of this technique to produce high quality thin 
films for many applications [I]. Those models adopt simple rules for the aggregation, 
diffusion, and desorption processes, consequently allowing the study of morphological 
properties of large samples Emu]. In the simplest cases, they assume limited mobility 
(LM) of adatoms; some examples are the models of Wolf and Villain and of Das 
Sarma and Tamborenea [^, in which short-range surface diffusion and permanent 
aggregation take place immediatly after adsorption. More realistic models consider 
thermally activated microscopic processes and are frequently called colletive diffusion 
models. The most prominent example is the Clarke-Vvedensky (CV) model [7], in 
which the adatom diffusion coefficients have Arrhenius forms, with energy barriers 
depending on the local number of lattice neighbors. An important difference from LM 
models is that the CV model obeys detailed balance conditions, thus it may also be a 
reliable description of the film dynamics without deposition Eli. 

In the basic formulation of the CV model, the dynamics may be represented 
by temperature-like parameters R and e, respectively representing the diffusion-to- 
deposition ratio of isolated atoms in terraces and the detachment probability at 
step edges i i. For the description of specific MBE processes, at submonolayer 
or multilayer regime, the models usually include additional energy barriers for 
diffusion across step edges (upward and downward movements) and additional adatom 
interactions; a thorough review of homoepitaxy applications is presented in Ref. [4] 
and recent advance is discussed in Refs. [g [ini nn [12]. The formation of patterns, 
growth and coarsening of mounds were some of the features that attracted much 
interest. 

However, a small number of works analyzed the dynamic scaling of surface 
roughness in CV-type models. This analysis provides a set of scaling exponents 
connecting the model to stochastic growth equations [2] , which helps to distinguish the 
essential physico-chemical mechanisms of film growth. Initial works on the basic CV 
model suggested temperature-dependent exponents and anomalous scaling of surface 
roughness [HlUlIilllS]. Subsequently, renormalization studies [nnnn suggested 
that it belongs to the class of the Villain-Lai-Das Sarma (VLDS) growth equation 
[201 [n]. In 2 -|- 1 dimensions, VLDS scaling was numerically confirmed for e = 0 and 
a broad range of R in Ref. [22] and for some values of R and e > 0 in Ref. [23]. 

The aim of this paper is to perform a systematic investigation of dynamic scaling 
in the basic CV model (without additional barriers across edges), with particular 
attention on the role of the temperature-like parameters. First, we study the scaling 
of global and local surface roughness and confirm that it belongs to the VLDS class for 
a broad range of values of R and e. Deviations from this scaling are shown to appear 
only when the surfaces are very flat. Second, we will show evidence of asymptotic 
normal scaling, but with an apparent anomaly for short times and small R similar 
to other VLDS models [ 22 ] • Finally, we will show that the roughness scales with R 
in a form similar to the irreversible aggregation model, with a weak dependence on 
the step detachment rate e. Although this model without barriers at step edges is of 
limited applicability to real solid films, these results may help the analysis of dynamic 
scaling in extended versions of the CV model, particularly if crossover features have 
to be analyzed. 

The rest of this paper is organized as follows. In Sec. [21 we present the model 
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and the related growth equation. In Sec. [31 we discuss the dynamic scaling of the 
global surface roughness. In Sec. HI the scaling of local surface roughness is analyzed. 
In Sec. ini we present our conclusions. 

2. Basic definitions and concepts 

2.1. Model and simulations 

The CV model is defined in a simple cubic lattice, with an initially flat substrate at 
z = 0. Deposition occurs with a flux of F atoms per site per unit time, in the z 
direction towards the substrate. Each adatom occupies one lattice site, whose side is 
taken as the unit length. We impose the solid-on-solid condition (i. e. overhangs are 
not allowed), thus only adatoms at the top of each substrate column are mobile. 

The hopping rate of an adatom with no lateral neighbor is 

Do = i'oex.p{-Es/kBT) ( 1 ) 

where vq is a frequency. Eg is an activation energy, and T is the temperature. The 
adatom step occurs in a randomly chosen substrate direction (Ex, Ey), towards the 
top of a NN column. If an adatom has n lateral neighbors, its hopping rate is 

D = Do€^ , e = exp(-Et/kBT), (2) 

where Ei, is a bond energy. 

An important parameter of the model is the diffusion-to-deposition ratio 

i? = ^ = y exp (-Eg/kBT). (3) 

It is usually interpreted as the number of steps of an adatom in a terrace before it 
is buried by the next atomic layer. However, it is highly probable that an adatom 
meets a lateral neighbor before being buried, which restricts that interpretation |22) . 
In the original CV model, vq = 2kBT/h, where h is the Planck’s constant, as predicted 
by transition state theory [7]. Some authors adopted that form [25j . but it is more 
frequent that a constant value t'o is considered in simulation and analytical 

works [4]. Here we will follow this trend and consider a fixed ratio vo/F = 10^^. 

Our simulations will be limited to deposition of 10'* monolayers. It is a typical 
value for thin films, corresponding to thicknesses of order 2 — 3^m, possibly more for 
molecular materials. The simulation time t will be expressed in number of deposited 
layers. The substrate size is L = 1024, which is large enough to avoid finite-size effects 
in the chosen deposition time. The values of Eg and E), are determined by material 
properties. Thus, since R = the parameters R and e simulataneously vary 

with the temperature for a given material. However, here we are interested in exploring 
a variety of physico-chemical conditions, which include different values of activation 
energies. For this reason, R and e will be taken as the independent parameters of the 
model. 

We will perform simulations in the range 10 < i? < 10*. Larger values are 
expected in many MBE processes, but we will show that they produce very smooth 
surfaces up to the maximal simulated thicknesses. We will also restrict our study 
to e < 0.1, since larger values of this parameter would represent a solid close to the 
melting point. 
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2.2. Dynamic scaling and universality classes 

The main quantity to characterize the film surface is the local roughness w (r, t) in 
boxes of size r at time t. For calculating this quantity, a square box of lateral size 
r glides along the film surface and, at each position, the root-mean-square (rms) 
height fluctuation of columns inside the box is calculated. The average among all box 
positions and among different configurations of the film at time t is the local roughness. 
The global roughness W (t) is measured in the full system size L, i. e. W (t) = w {L, t). 
In this work, very large substrates are considered, thus L has negligible effect on W. 

In systems with normal roughening (in opposition to anomalous roughening m), 
the expected scaling of the local roughness in large substrates is 

w(r,t) = (4) 

where a and z are the roughness and dynamic exponents, respectively, and / is a 
scaling function. For x = rjt^!^ <C 1 (small box sizes), g(x) is constant; for x ^ 1 
(large box sizes), the local roughness converges to the global one, W {t), which scales 
as 

W - t^, (5) 

where P = otjz is the growth exponent. 

When growth is dominated by surface diffusion, it is expected to be described by 
a fourth order stochastic equation in the hydrodynamic limit [ 2 ]: 

= + + T](f,t), ( 6 ) 

at 

where /i(r, f) is the height at position r and time t in a d-dimensional substrate, 
and A 4 are constants and 77 is a Gaussian, nonconservative noise (the contribution of 
the average external flux is omitted in Eq. El). The linear version (A 4 = 0) is the 
Mullins-Herring (MH) equation ^7}, while the nonlinear case is the VLDS equation 

MM- 

For the VLDS class in 2 -|- 1 dimensions, the best estimates of scaling exponents 
are given by the conserved restricted solid-on-solid models [35], and are very close to 
one-loop renormalization values [29]: a ~ 2/3, z ~ 10/3, and /3 ~ 1/5. 

3. Global roughness scaling 

In Fig. [T](a),we show the roughness evolution for R = 10^ and several values of e. The 
time scaling gives /3 0.20. As expected, the trend is that W decreases as e increases, 

since the detachment from steps of atoms with one or two lateral bonds helps them to 
move to positions with lower energy, forming more compact configurations. However, 
the quantitative effect of e on the roughness is small. 

In Fig. [1] (b), we show the roughness evolution for several values of i?, with 
e = 0.05. In order to estimate the exponent P, fits of the data are done only for 
VF > 2, because smaller roughness corresponds to very smooth surfaces. Those fits 
give 0.18 < P < 0.22, in good agreement with the VLDS exponent. For R = 10^, the 
roughness is very small at all times simulated here, thus deviations appear. 

As shown in Ref. [33], for e = 0, IT exceeds 2 units only at t ^ 10®; for e > 0, this 
occurs for longer times. Thus, for the values R > 10^ typical of most MBE processes, 
the film surfaces grown with the basic CV model are very flat. The formation of 
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Figure 1. (a) Global roughness as a function of time for R = 10^, with e = 0 (blue 
squares), 0.01 (red circles), 0.05 (black hexagons), and 0.1 (violet pentagons). The 
dashed line has slope 0.20. (b) Global roughness as a function of time for e = 0.05, 
with R = 10 (blue asterisks), 10^ (black hexagons), 10^ (orange stars), and 10^ 
(green squares). The dotted/dashed/full line has slope 0.22/0.20/0.18. 


patterns (e. g. mounds) observed in many simulations is possible only with energy 
barriers at step edges. 

Now we analyze the combined effects of parameters R and e on the roughness 
scaling. In the case e = 0, Ref. [22] showed that 




fO.2 




(7) 


in the growth regime. This was derived from a Family-Vicsek relation |30) that 
proposed the correlation length as ^ ^ (Rt)^^^ [due to the subdiffusive propagation of 
correlations; R (x D from Eq. (I3|)] and the saturation roughness as Ws ^ (due 
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to the formation of plateaus of width | 22 ) . 

Growth with e > 0 allows the detachment of adatoms from steps. It helps filling 
narrow surface valleys, which reduces the roughness. However, detachment may occur 
at upward and downward edges of the plateaus, thus is has a small contribution to 
their size. This explains why e > 0 does not lead to drastic changes in the roughness. 
Consequently, the aggregation of free adatoms at the step edges is still the main 
mechanism to determine the size of plateaus, and the dependence oi W on R is 
expected to be the same of the model with irreversible step aggregation [Eq. ([7|)] . 

The renormalization study of the CV model helps to infer the dependence of the 
roughness on the detachment probability. Ref. m derived the coefficients of the 
corresponding Langevin equation (with VLDS form) as a function of D and 7 = 1 — e. 
Those coefficients are products of D or by factors of the form {A + Bef' ^ where A 
and B are constants that depend on the coefficients of regularization of step functions 
(which are sensitive to discrete model rules) and c is an integer. 
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Figure 2. Global roughness as a function of scaled time according to Eq. m, 
with a = 0.025. The dashed line has slope 0.2. 


Following this reasoning, we propose a scaling relation for the CV roughness as 


W 


t 

i?3/2(e + a)®J ’ 


( 8 ) 


where 3 is a scaling function and a and 9 are constants. We did not restrict the scaling 
to the power-law regime [Eqs. ([5]) or ([ 7 ])] because this form may include the initial 
roughening (small W). 

Fig. [5]shows the roughness for several values of R and e with time scaled according 
to Eq. dS]). The excellent data collapse with 9=1 and a = 0.025 confirms the proposed 
scaling relation. For e < a, this relation clearly shows the weak dependence of W on 
that parameter; for e < 10“3 (£’f, > YksT), the effect of the detachment rate is 
negligible. The scaling of Eq. dH) also includes the case e = 0, in which detailed 
balance fails. 
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Figure 3. Local roughness as a function of box size for (a) i? = 10 , (b) R = 10^, 
(c) R = 10®, and (d) R = 10“^. Dashed (black) and full (red) lines indicate 
e = 0.01 and 0.1, respectively. For each set of parameters, from bottom to top, 
deposition times are i = 10, 10^, 5 X 10^, and 10®. For R < 10® and e = 0.01, 
t = 10^ is also included . 


The effects of R and e on the roughness are very different from their effects on 
island size in the submonolayer regime, in which scaling variables of the form Re'^ 
appear in the crossovers (with rational m > 1) [31j . Previous works on the basic 
CV model suggested deviations from VLDS scaling which can be addressed in the 
light of our results. The results of Figs. [T^ and [Di are similar to those shown in 
Ref. [TS], which worked in the range 10^ < R < 10^. That work suggested the 
possibility of non-universal exponents due to the small slopes of the fits of large R 
data, but, as discussed above, roughness scaling is not expected to appear in these 
conditions. Temperature-dependent exponents j3 were also suggested in Ref. |16j . 
which simulated the model with R < 10^ up to « 10^ monolayers. A possible crossover 
from uncorrelated deposition [ 2 ] (/3 = 0.5) to Edwards-Wilkinson scaling [ 32 ] (/? = 0) 
was suggested. However, the former was obtained for R ^ 1, in which most atoms 
cannot execute a single step, and the latter was obtained for very small values of the 
roughness, before the scaling regime. 

4. Local roughness scaling 

Figs. |3^-d show the local roughness for some values of R and e. The results for 
small R (100 or less) and small thicknesses (1000 monolayers or less) show features of 
anomalous roughening, i. e. the logic x logr curves for different times are split for 
small r. This was already observed for e = 0 in Ref. [33] . 
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Figure 4. (a), (b): time evolution of the local roughness for fixed box size, r = 5, 
respectively for = 10 and R = 10^; detachment parameters are e = 0 (blue full 
squares), e = 0.01 (black open squares), and e = 0.1 (red half squares); solid lines 
are linear fits with the effective anomaly exponents shown in the plot, (c) Local 
roughness at r = 5 as a function of for R = 10 and e = 0.01, using y = 0.09. 


For R > 10^, the curves for short times are also split, but coincide after « 100 
monolayers. The values of the local roughness are also very small, thus they would be 
hard to be distinguished in an experiment. Thus, the anomalous features disappear 
for typical MBE temperatures (i? > 10®); in this case, the presence of anomalous 
scaling in experimental data is a clear indication of the presence of energy barriers at 
step edges or other mechanisms that prevent surface smoothing by diffusion. 

Figs, m (a) and (b) show the time evolution of the roughness for fixed box size, 
r = 5, in growth with small R and thicknesses between 10 and 1000. Linear fits 
of those plots give effective anomaly exponents in the range 0.08 < k < 0.23. This 
procedure parallels the ones used in experimental works. It confirms that thin films 
grown at very low temperatures may show anomalous scaling features even in the 
absence of step edge barriers. 

However, this anomaly is only apparent, similarly to what occurs in other models 
in the VLDS class [24]. It is related to slowly vanishing (instead of increasing) terms 
in the scaling of the local slopes or of the small box local roughness, as 

w{ro,t) ^ A + (9) 

In Fig. m (c), we show w (5,t) versus for i? = 10 and e = 0.01, using y = 0.09. 
The good linear fit confirms the asymptotically normal scaling, but the small value of 
y gives large corrections even at long times [24] . 

The local roughness scaling does not provide reliable estimates of exponent a. 
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For small R, the scaling region of the logic x logr plot is very short (Fig. [3^) and the 
approximate slope is much smaller than a ~ 0.67 [28) . For large R, a longer scaling 
region appears, but the slope is also very small. On the other hand, the dynamic 
exponent z can be estimated by the procedure proposed in Ref. [33] . The first step is 
to calculate a characteristic length which is proportional to the correlation length 
at a given time t. This is obtained by defining rc as 

w {rc,t) = kW (t), (10) 

where W (t) is the global width and k is a constant. From Eqs. ([5]) and (|4]), it is 
expected that 

rc'-t^^^. ( 11 ) 

Here, we consider k = 0.7 for calculating Tc- 



1.0 1.5 2.0 2.5 

iogio(t) 


Figure 5. Time evolution of Vc for R=10 (triangles) and R=100 (squares). From 
bottom to top, detachment parameters are e = 0, 0.01, 0.1 for R=10, and e = 0, 
0.01 for R=100. The dashed lines have slopes 0.3. 


Fig. [S] shows the time evolution of Tc for some values of R and e; they are 
restricted to i? < 100 because the roughness was very small for larger R. Linear Ets 
give slopes between 0.30 and 0.32, which are in excellent agreement with the VLDS 
value 1/z ~ 0.30. It confirms that estimates of dynamic exponents from the local 
roughness are more reliable than those of roughness exponents [34j . 

The apparent anomaly for small R does not affect the estimates of z and /3, nor 
the global roughness scaling shown in Sec. (3] Due to these weak corrections, we 
propose an scaling relation for the local roughness considering: i) the general form of 
Family-Vicsek relation (|3]); ii) the same temperature-dependent variable (e-|- a) 
of the relation ([5]); iii) the scaling of the correlation length consistent with ^ ~ (Rt)^^^ 
for the case e = 0 [33] . This leads to 


t 

/3 

r 

_R^D (g -I- a)^ 

9 

(i?(e -1- 


w (r, t) = 


( 12 ) 
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where 5 is a scaling function. Fig. | 6 ] shows w/ \t/ (e + a))]^ as a function of 


r! i? (e + a)^^^t , for several values of i?, e, and t (data with very small roughness 


-|l/z 


is excluded, typically for short times or large R). The data collapse is also good, 
showing that Eq. (IT^ contains the leading temperature-dependent terms of the local 
roughness scaling. 



r/[R(E+a)“’t]''" 


Figure 6. Scaled local roughness as a function of scaled box size. 


5. Conclusion 

We performed simulations of the basic CV model (without extra energy barriers at 
edges) in 2 -f 1 dimensions for several values of the diffusion-to-deposition ratio R and 
detachment probability e in order to study the scaling of global and local roughness. 
Relatively short times consistent with thin film growth were considered. 

The exponents jd (z) calculated from the global (local) roughness confirm VLDS 
scaling with weak corrections. The scaling variable R^/"^ {e -\- a), with a = 0.025, 
represents the temperature effects in the dynamic scaling relations, with excellent 
accuracy for the global roughness (time scaling) and small corrections for the local 
roughness (time and box size scaling). This shows that R is the most important 
parameter to determine the surface morphology, with much smaller effects of e. This 
result for thin films is very different from that in submonolayer growth, in which the 
scaling of island size and related quantities combines rational powers of R and e [31j . 
On the other hand, the present scaling variables involving R and e are consistent with 
results of renormalization studies [T5] . 

The local roughness for R < 10^ shows evidence of anomalous scaling in the range 
of thicknesses considered here. However, asymptotic normal scaling is observed, with 
huge corrections expected even at long times. This is consistent with other VLDS 
models. 

Recent works on models with energy barriers for hopping across steps determined 
exponents /3 « 0.31 and 1/z ~ 0.22 [23l |35] for a range of temperatures; these 
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values are quite different from the VLDS exponents. The present model is of limited 
applicability to real solid films because they usually show this type of additional 
energy barriers. However, the methods proposed here may be useful for studying 
those extended models, for instance incorporating the effect of other temperature- 
dependent variables in the dynamic scaling, or investigating the question of anomalous 
versus normal roughening. 
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